DUDLEY KNOX LIBRARY 
NAVAL PC :J;TGRADUATE SCHOOL 
MONTERE ' 939434101 



Approved for public release; distribution is unlimited. 

Broadband Modal Beamforming 
of Acoustic Tomography Signals 
Acquired by a Vertical Array 

by 

Glenn A. Omans II 
Lieutenant, United States Navy 
B.S., Virginia Military Institute 

Submitted in partial fulfillment 
of the requirements for the degree of 

MASTER OF SCIENCE IN APPLIED SCIENCE 

from the 



NAVAL POSTGRADUATE SCHOOL 

Sept^^er 1992 yq 






Q^^i'15 

i.f 



ABSTRACT 

The objective of this thesis is to develop a technique 
and associated algorithms to extract the arrival time of modal 
energy, using a vertical array, from broadband signals. Modal 
energy arrival time is important to shallow water acoustic 
tomography because low angle rays, which contain the majority 
of acoustic energy, are often not resolvable. Tilt 
compensation is included in the beamforming algorithm to 
provide a virtual vertical array. A broadband modal filtering 
technique is accomplished through weighting the frequency 
components of phase encoded tomographic signals by the 
spectrum of the mode shapes. A methodology of phase decoding 
after beamforming was adopted to minimize processing. Initial 
development and prototyping was done using a parabolic 
equation model. Further testing was accomplished on real data 
taken from the Barents Sea Polar Front Experiment, August 
1992. Results show consistency over a number of transmitted 
pulses. Mode energy travel time measurement is simplified due 
to the distinct arrival structure of beamformed signals. 
Based on these results, the modal beamforming algorithm should 
be a useful tool for acoustic tomograpy. 
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I. INTRODUCTION 



A. BACKGROUND 

Tomography signal processing requires the precise 
measurement of acoustic travel time, which is the integral 
along the raypath of inverse sound speed. Travel time of 
acoustic signals has been well determined to be a function of 
temperature, salinity, and pressure. Once time perturbations 
for multiple arrival paths are known, ocean fluctuations can 
be determined using mathematical inverse techniques [Ref. 1] . 
In addition to measuring time perturbations based on ray 
arrival times, it is possible to measure perturbations based 
on mode arrivals. There are two approaches to tomography 
using modes: for continuous wave signals, the basic 

measurement is the modal phase difference; for broad band 
signals, the measurement used is modal travel time [Ref. 2] . 
Since coded signals (maximal -length sequence phase encoded 
signals) are intrinsic broadband signals, modal travel time 
measurement is germane to the work done in this paper. Where 
ray theory is not applicable, or not convenient, it is 
possible to extend tomographic techniques using both rays and 
modes [Ref. 3] . 

The goal of my thesis is to develop a technique for 
broadband modal beamforming for acoustic tomography. The 
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scope entails broadband modal beamforming with a vertical 
array and decoding (m- sequence removal of phase encoded 
signals) . In the process of my research several software 
programs were developed. Most notable is the broadband modal 
beamformer, which is an adaptation of the continuous wave 
beamformer developed by Crocker [Ref. 4]. Specific goals of 
research include: 

• Develop a broad spectrum modal filtering technique given 
the normal mode's eigen function and values. 

• Evaluate compatibility for phase sequence removal prior to 
and following beamforming. 

• Minimize processing time, where possible, and provide for 
a robust environment. 

• Evaluate performance using both synthetic and in situ 
data. 

• Provide a virtual array based on, time varying array tilt, 
array geometry and modal group speed. 

• Incorporate a supporting cast of programs, including a 
modal decomposition program for standardizing input to the 
beamformer and to simplify and enhance future signal 
processing. 

B. THESIS SUMMARY 

Chapter II, Barents Sea Experiment, is a summary of the 
Barents Sea Polar Front Experiment (BSPFEX) . This experiment 
is addressed because it provided real data to evaluate the 
processing algorithms. 

Chapter III, Acoustic Normal Mode Propagation Theory, is 
an introduction to normal modes. Modal acoustic propagation 
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is addressed. The chapter is concluded with a section linking 
acoustic ray and normal mode theory. 

Chapter IV, Signal Processing Concepts, presents the 
concepts and theory behind broadband modal beamforming. 
Subjects include beam steering, mode filtering and phase 
encoded signals. A plane wave beam steering approach is taken 
to compensate for array tilt. Following tilt correction, my 
modal filtering methodology is presented. Finally, a brief 
introduction to m- sequence signal decoding is given. Phase 
preservation after beamforming is addressed, supporting the 
concept of decoding after beamforming. 

Chapter V, Barents Sea Array, is a description of the 
array used during the Barents Sea Polar Front Experiment 
(BSPFEX) . Improvements made from previous vertical arrays are 
discussed. Statistics on hydrophone performance as well as 
the array modal beam patterns are addressed. 

Chapter VI, Evaluation of Barents Sea Data, is a 
presentation of results from the BSPFEX. Comparisons and 
conclusions are drawn concerning the array's performance. 

Chapter VII, Summary and Future Work, presents a summation 
of conclusions and identifies areas requiring improvement. 

The appendix contains the most significant algorithms I 
developed in the course of research. Algorithms included are: 
the broad band modal beamformer, modal decomposition program 
and the supporting modal extraction program which formats the 
input files for use in the beamformer. 
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II. BARENTS SEA EXPERIMENT 



Concepts developed in this thesis were applied to real 
data collected from the Barents Sea Polar Front Experiment 
(BSPFEX) . During the BSPFEX, I developed and tested my 
beamforming techniques and algorithms on synthetic data 
modeled using the Finite Element Parabolic Equation (PE) Model 
developed by Collins and Westwood [Ref. 5] . The synthetic 
data was modeled for the anticipated location of the array for 
the BSPFEX, Unfortunately, the array was deployed in a 
different location from the model limiting any comparisons 
drawn between the synthetic and BSPFEX data. 

The BSPFEX was conducted during August of 1992. The 
experiment comprised a joint effort between the Naval 
Postgraduate School (NPS) , Woods Hole Oceanographic 
Institution (WHOI) , and the Science Applications International 
Corporation (SAIC) . The principal investigators for the 
experiment are Professors Robert Bourke, Ching-Sang Chiu, and 
James H. Miller from NPS, Dr, James F. Lynch and Dr. A1 J. 
Plueddemann from WHOI, and Dr, Robin Muench from SAIC. The 
principal engineer for the vertical hydrophone array system 
was Mr. Keith Von der Heydt from WHOI. The objectives of the 
experiment as outlined by the Barents Sea Polar Front Group 
1992 are: 
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1. Provide a detailed physical description of the polar 
front . 

2. Enhance the understanding of dynamics of the front, 
including f rontogenesis and its influence on regional 
oceanographic processes. 

3. Assess the ability of acoustic tomography to define 
frontal and associated mesoscale features. 

4. Provide improved documentation of shallow water acoustic 
propagation in this region and the effect of the environment 
on acoustic ASW operations. 

[Ref. 6] 



This experiment is a milestone in acoustic tomography. It 
is the first time that a vertical array has been used in a 
shallow water environment to conduct acoustic tomography. A 
sixteen hydrophone telemetered vertical array, with ten meter 
spacing between hydrophones, was deployed during the 
experiment. The source and array mooring locations are listed 
in Table I. Table II describes the characteristics of the 
deployed broadband sources . 

Table I MOORING LOCATIONS FOR THE 224 AND 400 HZ SOURCES 
AND THE VERTICAL ARRAY. 



Mooring 


Latitude 


Longitude 


Depth 


SE(a) VLA 


74° 


19.1512'N 


23° 


33 .1438'E 


275.0 


SE(b) VLA 


74° 


19.1996'N 


23° 


32.2960'E 


275.0 


NE 224 Hz 


74° 


37.5535'N 


23° 


24.3755'E 


142.0 


NW 400 Hz 


74° 


32 .9152'N 


21° 


44.1043'E 


176.0 


SW 400 Hz 


74° 


04.7337'N 


22° 


00.4605'E 


380.0 
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Table II BSPFEX SOURCE CHARACTERISTICS. 



Mooring 


NE 


SW 


NW 


Freq (Hz) 


224 


400 


400 


SL (dB) 


183 


183 


183 


Duration (s) 


118.25 


132.86 


132.86 


BW (Hz) 


16 


100 


100 


Q (cyc/dig) 


14 


4 


4 


No. digits 


63 


511 


511 


Seq length 


3.9375 


5.11 


5.11 


No. sequences 


30 


24 +2 


24 +2 


Cycle (min) 2.5, 7.5 


... 0,10,20 ... 


5,15,25 ... 



Overall, the experiment was a sterling success. The WHOI- 
NPS vertical array was deployed on 12 August 1992. The array 
performed flawlessly for fourteen hours when a failure in the 
amplifier chip forced retrieval and repair of the array. 
After repair, the array was redeployed without the adjustable 
gain control making signal quantization from analogue to 
digital less optimal. Following recovery and redeployment, 
the array operated continuously for three days recording 
signals from both the 224 Hz and 400 Hz sources. In all, 
fifteen gigabytes of data were recorded. 

One major disappointment was the failure of the NW 400 Hz 
transceiver to transmit shortly after its deployment (see 
Figure 1) . No attempt was made to recover and repair the 
broken transceiver so as to not disturb the tomographic 
deployment schedule, and in the hopes that the transceiver was 
still able to receive acoustic data. Post experiment recovery 
of the broken transceiver revealed that it was still able to 
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and sources, from the BSPFEX. SW - 224 Hz transceiver; SE 
- 16 channel vertical array; NW/NE - 400 Hz transceivers 
[Ref. 7] . 

receive and recorded three days of acoustic data. With the NW 
transceiver still partially operational, three paths for 
tomography are available; two across the polar front and one 
along the front . 

In addition to tomographic pulse receptions, other 
acoustic obse^rvations were conducted. The ship dropped eight 
SUS charges to measure reverberation and bottom properties. 
Eighteen additional SUS charges were dropped from a P-3 
aircraft along the southwestern and northern tracks, to be 
used in evaluating propagation and bottom loss. One hour of 
continuous wave (CW) 224 Hz transmission was made for use in 
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modal phase and matched field tomography, and to measure 
temporal coherence for comparison with pulsed transmissions. 

Other sensors used during the deployment include, over 
three hundred Conductivity-Temperature-Depth (CTD) 
deployments, two of three current meters (one failure) , and 
acoustic doppler current profiler (ADCP) . Figure 2 denotes 
the fine structure of the polar front present, as measured 
from CTD deployments. Acoustic propagation was excellent 
during the experiment. The 224 Hz and 400 Hz sources were 
clearly audible at ranges over 50 kilometers, even without 
signal processing. With such clear receptions and low noise 
levels, the task of processing the data later will be much 
simplified. The above summary was based on the preliminary 
cruise report for the BSPFEX [Ref. 7] . 
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III. NORMAL MODE PROPAGATION THEORY 



A. ADIABATIC MODE THEORY 

Normal mode propagation theory is a physical model whose 
eigenfunctions, are solutions to the acoustic wave equation, 




at2 ' 



(3.1) 



in a waveguide [Ref. 8] . The objective of normal mode theory 
is to model the "local" pressure field as a linear combination 
of vertical standing waves. The term local refers to the 
range dependent sound speed structure and boundary conditions 
at a certain range. The following development applies for a 
slowly varying range dependent waveguide and is taken from 
Shang [Ref . 2] . 

The acoustic pressure field can be expressed in terms of 
local mode superposition. 



p(r, z) 



= E 

n=1 









(3.2) 



where '4'n(z,r) is the local solution or eigenfunction, and 
describes the distribution of modal energy as a function of 
depth. Fn(r) is the range dependent portion of the pressure 
field, and the term describes the geometric transmission 
losses caused by cylindrical spreading. 
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Subject to the boundary conditions and orthogonality, the 



local mode must satisfy. 



dz 








p (2) dz 






where k^(r,z) is the acoustic wave number, 3c„^(r) is the 
vertical wave number, or the eigen value, associated with the 
mode p(z) is the depth dependent density. The range 

dependent portion of the pressure field must satisfy. 



m 

dr^ 





~'E[2A^(r) 






, . dF„ 



(3.4) 



Ajnn Bnuj are coefficients describing the first and second 

order coupling effects in the sound channel: 



= / r—}^^iz,r ) — dz, 

J P^iz) dr 



Po(^) 



(3.5) 









d^^]f^(z, r) 
di^ 



dz. 



(3.6) 



If mode coupling is weak enough, then coupling can be 
neglected, and we have the adiabatic solution, 

Hr„(z|0)Ttf„(2|r) (3.7) 

v/^ 



p(r,z) = 

m 
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B. NOBMM. MODES AMD ACOUSTIC TOMOGRAPHY 



The travel time of each normal mode is dependent on the 
ocean structure. It is possible to estimate mode travel time 
from acoustic models, or it can be directly measured using 
mode filtering techniques. The difference between measured 
and estimated arrival times are the time perturbations for 
each mode. Broadband modal acoustic tomography seeks to 
exploit these time perturbations by applying inverse 
techniques to infer information about the acoustic 
environment . 

For a broadband source, the pressure field can be 
described as a function of a continuous wave field modified by 
the source signal spectrum. 



where s(o) is the signal spectrum of the source, and is the 
continuous wave field. The total pressure is the sum of the 
pressures for every mode. 




( 3 . 8 ) 



p(r, z, t) = z, t) , 



( 3 . 9 ) 



m 



where the modal pressure term is. 




( 3 . 10 ) 



12 



Taking the derivative of the phase with respect to frequency, 



t - 

"" dw 



r, 



(3.11) 



and it follows that the group speed must be, 



c ~ 



C?(0 



(3.12) 



It has been shown that there is a unique travel time 
associated with each mode for a given environment [Ref. 2] . 
Acoustic tomography can be accomplished using the time 
perturbations taken from the measured and estimated arrival 
times. The remaining task is to develop a technique to 
measure the modal arrival times, which is the subject of the 
following chapter on signal processing. 
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IV. SIGNAL PROCESSING CONCEPTS 



A. BROAD BAND MODAL BEAMFORMING 

Wave propagation modeling can be used to extract temporal 
and spatial characteristics for a propagating wave, and 
through inverse techniques to infer information about the 
acoustic medium. This is the essence of acoustic tomography. 
In recent years, ocean modeling for acoustic tomography has 
been dominated by acoustic ray theory. Ray theory is 
attractive because it is simple to model and requires minimal 
computer resources. Other modeling techniques simply were not 
practical until recent years. Ray theory is a high frequency 
approximation which fails to describe turning points, 
diffraction leakage and caustics. Further, many ray arrivals 
may not be resolvable, particularly in a shallow water 
environment . 

A vertical array provides an added dimension to 
tomography. Not only does it increase signal gain, it also 
allows mode and plane wave beamforming as well as enhancing 
resolvability of ray arrivals. The rest of this section will 
be devoted to the concepts required for modal beamforming. 
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1. Beam Steering 

Array tilt can be caused by a momber of forces such as 
waves, ocean currents, and external tensioning. Any tilt in 
the array is undesirable, and must be corrected. A virtual 
vertical array can be established using a plane wave beam 
steering methodology. 

There are essentially two approaches to plane wave 
beam steering. The first applies phase shift delays to the 
signal at the desired frequencies. The second uses true time 
delays . Phase delay beam steering is well suited for CW 
signals. However, this method is much less efficient for 
broad band signals because it requires phase shifting for 
every frequency of the broad band signal . 

Time delay beam steering is more suitable for use with 
broadband signals. Time shift delays are determined from 
array geometry and sound speed. For a line array, the time 
shift delay for the n^j element can be described as, 

, = , ( 4 . 1 ) 

" C 

where 6 ^ is the tilt in the direction of the source, and d is 
the element spacing [Ref. 9] . For normal mode beamforming, c 
is the group speed, equation (3.16), for each mode. Equation 
(4.1) is elementary in nature, and can be applied to any array 
or signal. Applying shift delays to discrete time signals 
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requires interpolation incurring additional processing [Ref. 
4] . 

2. Mode Filtering 

Plane wave beamforming does not account for the non- 
uniform energy distribution in the wave guide as described by 
normal modes. There are two types of modal beamforming, 
continuous wave (CW) and broad band (BB) . 

For CW modal beamforming, the task is simple. The 
time domain pressure field need only be weighted by the 
eigenfunction corresponding to each element depth. This 
method is not appropriate for broad band signals, because 
normal modes are a function of both depth and frequency. For 
BB beamforming, the frequency domain signal must be weighted 
by the corresponding frequency component of the eigenfunction. 
This requires that the mode shape be known at each frequency 
of the Discrete Fourier Transform (DFT) for the pressure field 
covering the desired bandwidth. 

For encoded signals with long duration sequence lengths, 
modal decomposition processing can be unsurmountable . The 
nominal modal decomposition time is thirteen minutes for each 
frequency. A 100 Hz bandwidth source with a sequence length 
of 5.11 seconds has 1022 frequency components. With current 
computer resources (SPARC 2) , it would take over nine days to 
process the eigen functions for each frequency. Fortunately, 
the mode shapes are continuous functions in depth and 
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frequency. Figure 3 shows how the mode amplitude, at a fixed 
depth, varies slowly with frequency. The mode shapes can be 
sampled as a function of frequency and then interpolated to 
achieve and desired frequency increment within the band. 




meter fixed depth. 



The next step is to weight the frequency domain 
pressure field by the mode shape at depths corresponding to 
the array elements. The following equation is the essence of 
modal filtering. 






= P ( 2) 









(4.3) 
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where is proportional to the intensity of mode m, and p is 
the density in the water column. The only limitation of 
equation (4.2) is that the Fourier Transforms are long enough 
to accommodate the length of the entire coded sequence. The 
pressure field frequency components are weighted by the 
frequency domain normal mode at the appropriate array element 
depths. The resultant is summed over depth, and an inverse 
Fourier Transform is applied to achieve the time domain 
beamformed signal. 

The array gain (AG) exhibited by broadband modal 
beamforming cannot be greater then conventional beamforming 
theory, 



AG = lOlog {N) , (4.3) 

[Ref. 10] . The theoretical array gain shown in equation (4.3) 
assumes the noise between elements is uncorrelated. Signal 
enhancement is strongly dependent on this correlation factor. 

B. MAXIMUM LENGTH PHASE ENCODED SIGNALS 

The objective of signal processing for broad spectrum 
acoustic tomography is the measurement of travel time for 
different raypaths and mode arrivals along a vertical cross 
section of the ocean. 

One method in which travel time can be measured is through 
the use of explosive and implosive devices. 
Unfortunately, explosive devices do not propagate 
repeatable acoustic signals. A better method that has 
been found employs the use of maximal -length sequences (m- 
sequences) or pseudo- random noise are well suited for this 
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application because of their deterministic nature, 
correlation properties, and simplicity [Ref. 11] . 

M- sequences are a sequence of binary bits having the property 

of a maximum possible period for an r-staged shift register. 

Shift registers are created from primitive polynomials. These 

binary shift registers are mapped to a series of ones and 

minus ones [Ref. 12]. The primary characteristic of the 

maximum length shift register is its autocorrelation 

properties. The self correlation of an m-sequence of length 

N produces a triangular function of short duration as shown at 

Figure 4 . 




Figure 4 Autocorrelation function for a m-sequence of 
length N. 



An m-sequence coded signal is formed using a simple 
harmonic source and applying a phase shift from the m-sequence 
register. A phase encoded signal has the form, 
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S(t) = ACOS + M( t) 4>) , 



(4.4) 



where A is the amplitude of the transmitted signal and M(t) is 
the value of the shift register for a specified digit 
duration. Q is an integer number of cycles (at carrier 
frequency) per digit and is used to determine the digit 
duration. The phase modulation angle <S> is chosen to optimize 
signal- to-noise performance; 

(|) = tan”^(/W) . (4.5) 

The one major drawback to using coded signals is that the 
standard correlation requires multiplications. Decoding 
algorithms can over come this processing issue by using the 
Fast Hadamard Transfoiw (FHT) . The FHT is analogous to a FFT 
and reduces processing to NlogN multiplications. 

Decoding signals is a simple process. The first step is 
to filter the input signal using a band pass filter covering 
the frequencies containing the coded signal. After filtering 
the signal is demodulated to remove the carrier frequency. 
The demodulated signal is then correlated for different shift 
versions using the FHT to reduce processing. Figure 5 is a 
block diagram showing the Birdsall -Metzger (BM) decoding 
process used in [Ref. 11] . 
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cos (2 Tt n/f 5 ) -*( X 



(7)«- sin(2uf ) 




Figure 5 Block diagram of m- sequence removal algorithm 
[Ref. 11] 
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C. SIGNAL PROCESSING SUMMARY 



The first step in modal beamforming is the correction of 
tilt to provide a virtual vertical array. This is 
accomplished using plane wave beamforming techniques. True 
time shift delays are used vice phase shift beamforming to 
reduce processing. After tilt correction, mode filtering is 
accomplished in the frequency domain. Modal beamforming 
constitutes filtering the frequency domain signals by the mode 
shapes, and then summing across the array elements. Post 
beamforming, the signals are decoded using the BM phase 
decoder algorithm. 

The phase decoding process can be accomplished prior to 
beamforming, but this would significantly increase signal 
processing as a multiple of the array elements and is not 
recommended. Further, the beamforming algorithm developed in 
this thesis does not support this procedure. If at some point 
it becomes desirable to process signals in this manner, the 
beamformer program would have to be altered to except complex 
data and to perform mode filtering about zero frequency offset 
signals . 
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V. BARENTS SEA ARRAY 



A. ARRAY COMPONENTS AND DESIGN 

The WHOI-NPS array contains sixteen hydrophone elements at 
ten meter intervals. The telemetered vertical array was 
moored to the bottom using a 1000 lb. cast anchor. The 
distance from the anchor to the bottom hydrophone is 
approximately 3 meters. A forty- eight inch bouyant syntactic 
foam sphere suspends the array 166 meters above its mooring. 
One meter below the sphere are the dual Benthos Interrogators. 
The interrogators were used to acoustically determine the 
exact location of the tethered array and to measure array 
tilt. Two inclinometers, one located two meters beneath the 
interrogator and the second near the center of the array, 
provided an additional means of determining array tilt. 
Figure 6 is an illustration of the array design. 

Nominal hydrophone sensitivity was -160 dB re 1 v/jwPa with 
a low frequency roll-off at 50 Hz. Data was preprocessed with 
an 8 -pole Butterworth low pass filter at 500 Hz cutoff 
frequency. A sampling rate of 1600 Hz was selected to achieve 
twice Nyquist frequency of the highest frequency source 400 
Hz. A sampling frequency four times carrier frequency is 
required by the m- sequence removal algorithm to simplify 



23 




Figure 6 Diagram of the WHOI-NPS vertical array used during 
the Barents Sea Polar Front Experiment [Ref. 7]. 



signal processing. Data from the 224 Hz source had to be 
resampled to 896 Hz sampling frequency for decoding. A single 
data bus was used to sample all sixteen array channels and 
created a 30 microsecond skew between elements. The 
beamformer corrects for the sampling skew based on the 
principles established for tilt compensation. 

B. BARENTS SEA POLAR FRONT EXPERIMENT - MODAL BEAMFORMER 
1. Equipment Performance 

The WHOI-NPS array was deployed in the Barents Sea on 
12 August 1992. The array performed well with one minor 
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mishap when the amplifier chip failed. The amplifier chip 
failure effected adjustable gain control and degraded analogue 
to digital quantization for much of the data collected. All 
array elements operated without failure, and the array 
remained operational throughout its deployment. 

A new mooring and tether design gave the array 
increased stability. The mean measured array tilt was three 
tenths of a degree. This is meaningful considering that 
previous vertical arrays suffered tilt in excess of ten 
degrees thereby degrading system performance. 

2 . Modal Decomposition 

Mode filtering performance is highly dependent on 
array geometry and the normal mode shapes. Modal 
decomposition requires knowledge of upper and lower boundary 
conditions, sound velocity, and density. The algorithm used 
for modal decomposition is located at the appendix of this 
thesis. [Ref. 13] The decomposition algorithm employs a 
finite differencing routine to approximate the initial eigen 
value problem. Implicit to the routine is the assumptions of 
a pressure release surface and a perfectly rigid lower 
boundary. The bottom is modeled using constant sound speed 
and density characteristics. There is no limitation to the 
quantity or complexity of bottoms except for processing time, 
which increases as the square of the number of points in the 
profile . 
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Sound and density profiles used to calculate 
eigenvalues and shapes for the Barents Sea were determined 
from CTD data, and estimations of bottom sound speed and 
density. All processing is based on a single bottom model 
with a constant sound speed of 1667 m/s, and density of 1.83 
g/cm^, taken from Clay and Medwin [Ref. 14] . 

Figure 6 is the Barents Sea sound velocity profile at 
the location of the array. The sound velocity indicates the 
existence of a well defined mixed layer followed by a strong 
negative sound speed gradient extending to the bottom. Due to 
the negative sound speed gradient, the majority of acoustic 
energy is expected to be trapped well below the surface with 
significant bottom interaction. 

3 . Normal Modes at the Vertical Array 

Figures 7 to 12 are the mode shapes for the Barents 
Sea array at the 224 Hz source center frequency. Modes one 
through four have turning depths well beneath the surface and 
are expected to contain the majority of acoustic energy based 
on the strong negative sound speed gradient. Modes five 
through ten also turn beneath the surface, but should have 
less beamformed energy than the first four modes because there 
is significant modal energy outside the array's aperture. 
Modes ten through twenty interact with the surface. These 
modes should have much less energy due to surface scattering 
effects . 
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Figure 7 Sound speed profile at the receiving array. 
Bottom sediment starts at 275 meters modeled by constant p 
1.83 g/cm^, and sound speed 1577 m/s; The bottom is assximed 
rigid below 420 meters . 

4. Modal Beam Pattern 

There are many measurements of performance for 
classical arrays. Plane wave beamforming is a well understood 
process, and performance of a plane wave beamformer can be 
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Figure 8 Normal modes 1 through 4 . 
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Figure 9 Normal modes 5 through 8 . 
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Figure 10 Normal modes 9 through 12 . 
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Ficpure 11 Normal modes 13 through 16. 
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Figure 12 Normal modes 17 through 20. 
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easily determined through a multitude of available measures 
including: theoretical array gain, estimated array gain, and 
the steering beam pattern. Array gain experienced by modal 
beamforming cannot exceed the theoretical gain for a classical 
beamformer. However, there are additional issues to address 
concerning tomographic resolvability for mode arrivals. 

One such factor, is the cross talk, or the recognition 
differential between modes. To quantify this affect, 
synthetic signals modeled from individual modes were 
beamformed. The modal signals are created using the frequency 
domain eigen solutions at each hydrophone depth. To duplicate 
the coded signals frequency structure, a Blackman Window is 
applied to the modal signal's frequency components. The 
square of the beamformed signal is proportional intensity of 
the mode. The peak values of the beamformed modes are squared 
and normalizing by the corresponding modal signal, forming the 
modal recognition kernel. The transpose of the recognition 
kernel gives the modal beam pattern. 

Figures 13 to 18 are the modal beam patterns specific 
to the Barents Sea array geometry and ocean structure. The 
following is a summary of my conclusions and conjectures on 
the Barents Sea array modal beam pattern; 

• Modes one through three are well resolved. 

• Modes four and above are not spatially well resolved due 
to the size of the array's aperture. 
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• Above mode three resolvability must be handled on a case 
to case bases, taking into account the arrival structure 
and relative intensity of each mode. 

• Assuming there is no appreciable modal energy above mode 
twenty, spatial aliasing effects are non-existent for this 
array . 

• Above mode four, there is a recurring beam pattern which 
indicates that the closest modes are the hardest for the 
beamformer to differentiate. Mode four beam pattern shows 
the transition from well resolved to less resolvable 
arrivals . 

• Increasing array length to 50 meters below the surface 
would increase the coverage of the array's aperture and 
guarantee that at least the first six modes are spatially 
resolved. This can be accomplished through increasing the 
number of hydrophones, or by expanding element spacing, 
subject to spatial aliasing. 
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Figure 13 Recognition beam pattern for modes one through 
four . 



35 



mode # mode # 



Mo d e 5 - Beam Pallcrn ,, Mode 6 - Beam Pallern 





Csl 



LO 






9P 



ap 




Figure 14 Recognition neam pattern for modes five through 
eight . 
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Figure 15 Recognition beam pattern for modes nine through 
twelve . 
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Figure 16 Recognition beam pattern for modes thirteen 
through sixteen. 



38 




ap 



ap 




ap ap 

Figure 17 Recognition beam pattern for modes seventeen 
through twenty. 
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VI. EVALUATION OF BARENTS SEA DATA 



A. DECODING OF RAW ACOUSTIC DATA 

The participants of the Barents Sea Polar Front Experiment 
performed some near- real -time signal processing during the 
experiment. Some of the signal processing included decoding 
of the 224 Hz source. Figure 18 is the decoded signal taken 
from the bottom hydrophone of the array. 

xlO-* Acoustic dme series after decoding 




Figure 18 Acoustic time series after BM decoding for the 
bottom hydrophone in the array [Ref. 7] . 

The BM sequence removal algorithm performed very well. 
Theoretical decoding gain is based on the ntjmber of digits in 
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the coded signal . The 224 Hz source signal has a theoretical 
decoding compression gain of 18 dB which compares favorably to 
initial estimates. Signal gain estimates where determined by, 



[Ref. 15] . 

Coherent averaging is a useful method for determining the 
acoustic arrival structure across the array. Figure 19 is the 
coherently averaged signal across the array, post m- sequence 



presence of at least two separate arrivals. This effect is 
most likely do to the time spreading of the signal between the 
propagation paths. The first arrival has the bulk of acoustic 
energy and a complex structure. Clearly this energy is formed 
by the bottom bounce propagation path. The second arrival 
lacks the energy and fine definition of the first caused by 
increased bottom and surface interaction. 

B. RESULTS FROM THE BROADBAND MODAL BEAMFORMER 

The results from broadband modal beamforming of the 
Barents Sea data are quite impressive. The mode arrival 
shapes are very distinctive, which simplifies arrival time 
estimation. All times for Figures 21-27 are relative to the 
record tape, and do not represent the travel from source 
emission. Figure 20 shows the mode one arrival structure for 




( 6 . 1 ) 



removal . The coherently averaged signal indicates the 
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After decoding and coherent average 




Figure 19 Coherent averaging after BM decoding for the 
WHOI-NPS array [Ref. 7] . 

the first four sequences of decoded pulses. Figures 21 
through 23 are the first five modes pertaining to the first 
four sequence arrivals of Figure 18. The first five modes are 
the strongest arrivals. This is expected because the majority 
of modal energy is within the array aperture. Mode three is 
the strongest arrival indicating the majority of acoustic 
energy is in the lower 125 meters of the water column. Mode 
two is particularly interesting due to its multiple arrivals. 
This pattern might be attributed to coupling caused by the 
sharp change in sound speed gradient at the front. The mode 
two arrival is both the sinallest in total energy and the 
earliest arrival of the first five beamformed modes. 
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Figure 20 Mode one, beamformed and decoded arrival 
structure . 

Modes six through ten, Figure 25, appear after the initial 
arrival structure. These inodes contain less energy and have 
later arrival times. A large amount of modal energy exists 
outside the array aperture and is at least a partial factor to 
the decrease in beamformed energy. 

Modes eleven through twenty. Figures 26 and 27, show a 
sharp drop off in energy levels which cannot be completely 
contributed to the array aperture. This sudden decrease of 
beamformed energy is most likely caused by high transmission 
loss caused by surface scattering. It is interesting to note 
that arrivals for modes sixteen and seventeen indicate 
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increased energy levels. This is a curious result that 
presently eludes explanation. 

Broadband modal beamforming of the BSPFEX data provided 
extremely clear arrival structure that could easily be 
mistaken for modeled data. Part of the reason that the 
results are this good is the extremely high signal to noise 
ratios and the low propagation loss. Clearly, the first ten 
modes are the most significant to acoustic tomography. Mode 
resolvability is dependent upon the mode's arrival time, the 
modal beam pattern and spatial aliasing. Resolvability issues 
have not yet been completely sorted out but may indicate that 
only the first few modes are reliable enough to use in the 
tomographic inverse . 
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21 First arrival sequence for modes one through five. 
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Figure 22 Second arrival sequence tor modes one through five. 
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^gure 23 Third arrival sequence for modes one through five. 
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Figure 24 Fourth arrival sequence for modes one through five. 
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Figure 25 Arrival structure for modes six through ten, 
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Figure 26 Arrival structure for modes eleven through fifteen 



cc 







VO 






c/: 

o 



ro 






Figure 27 Arrival structure after beamforming for modes 
fifteen through twenty. 
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VII. SUMMARY AND FUTURE WORK 



A. CONCLUSIONS AND SUMMARY 

Preliminary results indicate that broadband modal 
beamforming is a useful tool in shallow water acoustic 
tomography. Ray acoustic tomography in shallow water is 
complicated because the earliest ray arrivals, containing the 
bulk of the acoustic energy, are unresolved. Unresolved ray 
paths are well described by the lower mode arrivals of the 
beamformed array. Broadband modal beamforming produces 
distinct mode energy arrivals simplifying estimation of 
arrival times. 

Improvements made by WHOI and NPS to the vertical array 
design proved very successful . The Dual Benthos interrogators 
provided a more reliable source for measuring array geometry 
than previous methods. Array stability was vastly improved by 
its new moored design and should be used as a benchmark for 
the design of future arrays. The measured array tilt during 
the experiment was a nominal three tenths of a degree and 
minimized the correction required for array tilt. 

Based on recommendations from Crocker [Ref. 4] , a fast 
third order interpolation routine was implemented to reduce 
processing time. The new interpolator was derived using the 
Taylor Approximation for a third order polynomial [Ref. 16] . 
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This algorithm is approximately five times faster than 
previous interpolator based on Neville's algorithm 
[Ref. 17] . 

B. RECOMMENDATIONS 

The broadband modal beamformer requires each mode to be 
processed individually. Most of the time incurred from 
processing is due to the Fast Fourier Transform algorithm. 
Processing all modes simultaneously would reduce the number of 
Fast Fourier Transforms by the number of desired modes. 
Significant savings in processing can be realized by adopting 
this technique. Calculating the modal beam pattern was an 
additional time sink. Since each mode pattern must be 
processed for every mode, processing time goes as the square 
of the number of modes desired in the beam pattern. For 
twenty modes, the beamformer must be run four hundred times to 
get the beam pattern for just one sound velocity profile. If 
all the modes where processed simultaneously, the modal beam 
pattern could be an added option to the program reducing the 
processing required by the user. 

Decoding of the m- sequence coded signal occurs after 
beamforming. The decoding process could have been 

accomplished prior to beamforming but would have increased 
processing time by the number of elements in the array. The 
output from the m- sequence removal algorithm is complex 
numbers thereby doubling the amount of memory required to 
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store data. Despite this, it may still be desirable to decode 
the data first so that plane wave and modal beamforming can be 
accomplished using the same data base. To perform sequence 
removal prior to beamforming, the decoder algorithm must be 
adapted to handle multi-channel data, and an interpolation 
routine should be added to resample data prior m- sequence 
removal. The modal beamformer will also need to be updated to 
handle complex input and to perform mode filtering on 
demodulated signals. 

Design of the next vertical array should incorporate 
estimations of the modal beam pattern. Modeling of the 
acoustic environment will give an estimate of the modal beam 
pattern. Using this information, the number of array elements 
and geometry can be optimized to achieve the desired 
performance level . 
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APPENDIX 



A. BROADBAND MODAL BEAMFORMER ALGORITHM 

The broadband modal beamformer was designed for a robust 
environment. Program portability allows it to be used on any 
computer system supporting C or quick C software. The program 
is flexible enough to cope with most vertical line array 
designs. A large number of user parameters are incorporated 
in the computer architecture to enhance portability. The 
following is a summery of user parameters employed. 

1 . Program Parameters 



• UNIX VERSION: a compiler flag used to determine the 

operating system, either UNIX or ANSI. 

• ASCII: ASCII and BINARY options determine the beamformed 
output type. These options are inter-linked, and never 
can be activated at the same time. 

• SIGNAL: enables time series beamformed output. This 

parameter should be left "ON" . 

• LOWER_SENSOR: enables lower sensor tilt data if available. 

• VALIDATE: this option flags the program to output several 
validation files. If selected, the program will output 
files containing array geometry, steering delays, 
hydrophone weights at the lowest frequency in the band, 
and estimated mode speed. 

• ERROR_ESTIMATE : flags to produce the estimated error 

caused by interpolation. The error estimate is stored in 
a file specified by the user. 

• VERBOSE: if selected, this option will cause the program 
to provide a status report during operation. Information 
is sent to standard output. This option is useful for 
debugging and for short data sets. 

• ON/OFF: logical switches for program control. 
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• INTERPOLATE; specifies which interpolation routine to use. 
The available interpolators are fastpolyS, polint and 
ratint. FastpolyS is the fastest of the three 
interpolators. Other interpolators can be added by the 
user if desired. 

• ORDER: determines the order for interpolation. Polint and 
ratint are multi order algorithms, while fastpolyS can 
only perform third order interpolation. If the order of 
interpolation is not three and the fastpolyS interpolator 
is selected, the program will terminate on a standard 
error. 

• STEP: indicates the number of points on either side of a 
sequence for derivative estimates. 

• TINY: smallest number allowed for floating point 

operations . 

• PI : trigonometric constant . 

• RADIAN: degree to radian conversion constant. 

• OFFSET; correction term for the difference between tilt 
sensor depth and the depth of the first hydrophone. 

• DELTA_R: array element spacing in meters. 

• CTD_OFFSET: used to correct for offset caused by the CTD 
measurement. This offset is reflected in the eigen 
solutions . 

• SAMPLE_DELAY : corrects for sampling skew caused by 

sampling data from a single data bus (ms) . 

• SSP_LENGTH: dimensioning term based on the maximum number 
of points in a single eigen function. 

• EIGVAL_LENGTH: dimensioning term based on the maximum 

number of frequencies contain in the eigen value file. 

• LOOK DIRECTION: desired direction for beamforming in 

degrees true. 

• TILT_BUFFER: dimensioning factor based on the maximum 

number of tilt observations in the tilt input file. 

• BUFFER_TIME : dimensioning term. Must be an integer and 
greater then the FFT_TIME + 2 . 

• F_SAMPLE: the integer sampling frequency. 
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• CHANNELS : total number of array elements . 

• F_CARRIER; the desired carrier frequency to beamform. 
Note that this should match the center frequency of the 
eigen function and value input files. 

• FFT_TIME: amount of data, in seconds, allotted to each 
FFT. This should be the sequence length of the desired 
beamformed signal . 

• FFT_LENGTH: the radix two FFT size. The FFT_LENGTH must 
be larger or equal to FFT_TIME*F_SAMPLE . 

• SWAP: macro used by the FFT algorithm. 



2 . Beamformer Source Code 



j -k •k -k -k 

* PROGRAM: BEAMFORMER vsn 5.0 

* NAME : BBbeaxn . c 

* WRITTEN BY: Glenn A. Omans II / Steven Crocker 

k 

* PURPOSE: Broadband Modal Beamforming of a Vertical Array. 

* 

* LAST UPDATE: September 28, 1992 



* This program takes input from various data files and the user. It 

* outputs a data file. The inputs are a number of channels of digital 
acoustic data, and information regarding the physical characteristics 

* and geometry of the receiving array. Additionally, environmental data 

* in the form of normal mode eigenfunctions and eigenvalues at the 

* receiving array are required to operate this beamformer. The output 

* is a single channel of acoustic data. 



kkkk ^ 

#define UNIX VERSION /* either ANSI or UNIX */ 
#define ASCII ON /* select output mode ON or OFF */ 
#define BINARY OFF /* select output mode ON or OFF */ 
#define SIGNAL ON /* either ON or OFF */ 
#define LOWER_SENSOR OFF /* either ON or OFF */ 
#define VALIDATE OFF /* either ON or OFF */ 
#define ERROR_ESTIMATE OFF /* either ON or OFF */ 
ttdefine VERBOSE OFF /* either ON or OFF */ 
#define ON 1 /* logical "switch” */ 
#define OFF 0 /* logical "switch” */ 
#define INTERPOLATE fastpolyS /* either polint, fastpolyS or ratint */ 
#define ORDER 3 /* order of interpolator (odd)*/ 
#define STEP 1 /* number of steps for derivatives */ 
ttdefine TINY l.Oe-25 /* prevents division by zero */ 
#define PI 3.14159265359 /* for freq to omega conversions */ 
#define RADIAN 57.2957795131 /* for degree to radian conversions */ 
#define OFFSET 0.0 /★ dist btwn upper sensor and phone #1 */ 
#define DELTA_R 10.0 /* array element spacing */ 
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#def ine 
#def ine 
#def ine 
#define 
#define 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 



CTD_OFFSET 0.0 
SAMPLE_DELAY 0.003 
SSP_LENGTH 2500 
EIGVAL_LENGTH 230 
LOOK_DIRECTION 270.0 
TILT_BUFFER 120 
BUFFER__TIME 7 
F_SAMPLE 1600 
CHANNELS 16 
F_CARRIER 224.0 
FFT_TIME 3.9375 
FFT_LENGTH 8192 
SWAP (a,b) tempr= (a) ; (a) 



/* diff btwn ctd depth inc & 1st depth * 
/* array sampling delay (ms) ’ 

/* max number of pts in eigunfunction ’ 
/* max number of eigenvalues ’ 

/* direction from which signal arrives ^ 
/* max length of tilt data vectors * 
/* input buffer length in seconds (int)"* 
/* sampling frequency (int)^ 

/* number of channels processed * 

/* carrier frequency * 

/* time alloted to each fft * 

/* radix 2 <= (FFT_TIME-2) *F_SAMPLE * 
= (b) ; (b) =tempr 



#include<stdio .h> 
#include<malloc .h> 
#include<math . h> 



#if defined ( ANSI ) 

#include<f loat .h> 

#include<stdlib.h> 
int getinput (void) ; 
int putOutput (void) ; 

int processTilt (float **x, float **y, float **z, float *delayClock) ; 
int processModes (float **z, float **weight, float *ptrC) ; 
int dydx (float *x, float *y, float *ddx, int points); 

int fastpoly2 (float *xa, float *ya, int n, float x, float *y, float *dy) , 

int fastpoly3 (float *xa, float *ya, int n, float x, float *y, float *dy) , 

int polint (float *xa, float *ya, int n, float x, float *y, float *dy) ; 

int ratint (float *xa, float *ya, int n, float x, float *y, float *dy) ; 

int realft (float *data, int n, int isign) ; 

int f ourl (float *data, int nn, int isign) ; 

int window (float *data, int N) ; 

float *vector(int length); 

float **matrix(int row, int col) ; 

int free_matrix (float **m, int row) ; 

void ExitOnError (char error_txt [] ) ; 

#elif defined ( UNIX ) 
int vector 0, 
matrix ( ) , 
getinput () , 
putOutput ( ) , 
processTilt () , 
processModes () , 
dydx ( ) , 
f astpoly2 () , 
fastpoly3 () , 
polint 0 , 
ratint () , 
dumpSpectrumO , 
window ( ) , 
realft () , 
f ourl 0 , 
free_matrix 0 , 

ExitOnError 0 ; 

#endif 



/* Global Variable Declarations */ 
float **inSOUND, *outSOUND, *MeanSqError , Max=0.0, Min=1.0e25; 
float fmax, fmin; 

int LastDelay, Minute=l, f irstBuf f er=l ; 
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FILE *fpInSound, *fpOutSound, *fpOut Spectrum, *fpMSE; 



main ( ) 

{ 

char fileDelay[80] , f ileModes [80] , f ileArray [80] ; 

FILE *fpVl; 

int i, j, k, n, itime, ibuff, kount, *shift; 

float **x, **y, **z, *indx, *samples, arg, ans, err, *delay, 
*delta, **weight, *pwrSpectrum, Cgroup, *Xf, 

*sumXf, **buffSOUND, df, freq, *delayClock; 

double Time; 



/* Memory Allocation and Tilt Data Processing */ 
x= (float**) matrix (CHANNELS, TILT__BUFFER) ; 
y= (float**) matrix (CHANNELS, TILT_BUFFER) ; 
z= (float**) matrix (CHANNELS, TILT_BUFFER) ; 
delayClock= (float*) vector (TILT_BUFFER) ; 

processTilt (x, y, z, delayClock) ; 

/* printf ( ”\nDelay Clock Times\n\n") ; 

for (i=0; i<=LastDelay ; i++) printf ( ”%i\t%g\n” , i, delayClock [i] ) ; */ 

if (VALIDATE ==0N) 

printf (”\nEnter file name for array validation data: ”) ; 
if ( (scanf ( ”%s" , f ileArray) ) ==E0F) 

ExitOnError ( "Fatal error in scanfO"); 

fpVl=fopen (f ileArray, ”wt") ; 

if (fpVl==NULL) ExitOnError ( "Error opening validation file"); 

fprintf (fpVl, ”Channel\t\t X\t\t Y\t\t Z\n") ; 

for (i=l ; i<=LastDelay; i++) 

{ 

fprintf (fpVl, "MINUTE: %i\n" , i) ; 
for (j=l; j<=CHANNELS; j++) 

fprintf(fpVl, "%i\t %f\t %f\t 
%f \n" , j ,x [j] [i] ,y [ j] [i] , z [j] [i] ) ; 

} /* for */ 
f close (fpVl) ; 

} /* if */ 

/* Memory Allocation*/ 

weight= (float**) matrix (CHANNELS, FFT_LENGTH) ; 
indx= (float*) vector (ORDER+1) ; 
delay= (float*) vector (CHANNELS) ; 
delta= (float*) vector (CHANNELS) ; 

OUtS0UND= (float*) vector ( (BUFFER_TIME-1) *F_SAMPLE) ; 

Xf= (float*) vector ( (FFT_LENGTH+1) *2) ; 

BumXf = (float*) vector ( (FFT_LENGTH+1) *2) ; 

buffSOUND= (float**)matrix(CHANNELS, (BUFFER_TIME-1) *F_SAMPLE+1) ; 
inSOUND= (f loat**)matrix(CHANNELS, BUFFER_TIME*F_SAMPLE) ; 
if ( (shift= (int*)malloc( (CHANNELS+1) *sizeof (int) ) ) ==NULL) 

ExitOnError ("Memory allocation failure for shif t [] ."); 
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0 N ) 



for (i = l; i<=ORDER+l; i + + ) indx [i] = (float) i ; 

if (ERROR_ESTIMATE = = 

MeanSqError= (float*) vector ( (BUFFER_TIME-1) *F_SAMPLE) ; 

/* Mode Data Processing */ 
processModes ( 2 , weight, &Cgroup) / 



*/ 



/* Calculate delays, shifts, ect */ 
for (i=l; i<=CHANNELS; i++) 

{ 

delay [i] =x [i] [Minute] /Cgroup + /* Time delay */ 

(float) (i-1) *SAMPLE_DELAY/1000. ; /* Sampling delay */ 

shift [i] = (int) (delay [i] * (float) F_SAMPLE) ; /* # of samples */ 

/* fraction of 1 sample 

delta [i] =delay [i] * (float) F_SAMPLE- (float) shift [i] ; 

} /* for */ 



while (Minute<=LastDelay) 

{ 

if (VALIDATE==ON) 

{ 

if (f irstBuff er) 

{ 

printf ("\nEnter file name for delay validation: ") ; 
if ( (scanf C'%s",fileDelay) )==EOF) 

ExitOnError ("Fatal error in scanf () "); 

if ( (fpVl=fopen (f ileDelay, "wt")) == NULL) 

ExitOnError ( "Error opening validation file."); 

} /* if */ 
else 

^ if ( (fpVl=fopen(f ileDelay, "at")) == NULL) 

ExitOnError ("Error opening validation file."); 

} /* else */ 

fprintf (fpVl, "MINUTE: %i\n", Minute); 
fprintf (fpVl, 

"Channel\t delay\t\t int shift\t fraction of 1 shift\n") ; 

for (i=l;i<=CHANNELS;i++) 

fprintf (fpVl, "%i\t %e\t %i\t% e\n", 
i , delay [i] , shift [i] , delta [i] ) ; 

} /* for */ 
fclose(fpVl) ; 

if (f irstBuffer) 

printf ( "\nEnter file name for phone weights and group speed: 

") ; 

if ( (scanf ("%s" , f ileModes) ) ==EOF) 

ExitOnError ("Fatal error in scanf () "); 

if ( (fpVl=f open (f ileModes, "wt")) == NULL) 

ExitOnError ( "Error opening validation file.\n"); 
fprintf (fpVl, "Group speed for F_CARRIER is: 
%g\n\n\n", Cgroup) ; 
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} /* if firstBuffer */ 
else 
{ 

if ( (fpVl=fopen(fileModes, "at") ) == NULL) 

ExitOnError { "Error opening validation file.\n"); 
fprintf (fpVl, "\nHydrophone weights for minute %i\n" , Minute ) ; 
} /* else if not firstBuffer */ 

fprintf (fpVl, "Channel \tWeight\n") ; 

for (i=l;i<=CHANNELS;i++) 

fprintf (fpVl, "%i\t%e\n" , i, weight [i] [1] ) / 

f close (fpVl) ; 

} /* if VALIDATE */ 



Time =1; /* First second of time is 

chopped */ 

for (n=l; n<=60/ (BUFFER TIME-2); n++) /* for one Minute */ 

{ 

getinput () ; 



if (firstBuffer) /* Produce first output buffer */ 

if (VERBOSE) printf ( "Interpolation first Buffer\n"); 
for (i=l; i<=F SAMPLE; i++) 

{ 

OUtSOUND [i] =0.0; 

if (ERROR_ESTIMATE==ON) MeanSqError [i] =0 . 0 ; 

} /* for i */ 

for (i=F_SAMPLE+l; i<= (FFT_TIME+1) *F_SAMPLE ; i++) 

{ 

OUtSOUND [i] =0 . 0 ; 

if (ERROR_ESTIMATE==ON) MeanSqError [i] =0 . 0 ; 

Time+= (1/ (float) F_SAMPLE) ; 

if (Time >= delayCloclc [Minute] ) 

{ 

Minute++; 



delay */ 
Sampling delay 
samples */ 



/* Mode Data Processing */ 
processModes (z, weight, &Cgroup) ; 



/* Calculate delays, shifts, ect */ 
for (k=l; k<=CHANNELS; k++) 

{ 

delay [Jc] =x [k] [Minute] /Cgroup + 



*/ 



(float) (k-l) *SAMPLE_DELAY/1000 . ; 
shift [k] = (int) (delay [k] * (float) F_SAMPLE) ; 



fraction of 1 sample */ 



/* Time 
/* 

/* # of 
/* 



delta [k] =delay [k] * (float) F SAMPLE- (float) shift [k] ; 
} /* for */ 

} /* if (Need new shift delays) */ 
for (j=l; j<=CHANNELS; j++) 
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if (delay [j] >=0 . 0) 

{ 

arg= (float) (ORBER+l) /2 .0+ (l-delta[j] ) ; 
samples = &inSOUND[j] [i-shift [j] - (ORDER+1) /2] ; 

} /* if */ 

else if (delay [j] <0 .0) 

{ 

arg= (float) (ORDER+1) /2 . 0+delta [j] ; 

samples = &inSOUND [j] [i-shift [j] - (ORDER+1) /2+1] ; 

} /* else if */ 

INTERPOLATE (indx, samples, ORDER+1, arg, &ans, &err) ; 
buffSOUND[j] [i-F_SAMPLE] =ans; 

if (ERROR_ESTIMATE==ON) 

MeanSqError [i] =MeanSqError [i] +err*err; 

} /* for */ 

if (ERROR_ESTIMATE==ON) 

MeanSqError [i] =MeanSqError [i] / (float) CHANNELS ; 

} /* for */ 

} /* if */ 

else /* Produce subsequent output buffers */ 

{ 

if (VERBOSE) printf ( "Interpolating subsequent output 

buffers\n") ; 

for (i=F_SAMPLE+l; i<= (FFT_TIME+1) *F_SAMPLE ; i++) 

{ 

OUtSOUND [i-F_SAMPLE] =0.0; 

if (ERROR_ESTIMATE==ON) MeanSqError[i-F_SAMPLE]=0.0; 

Time+= (1/ (float) F_SAMPLE) ; /* time processed with first 

second chopped */ 

if (Time >= delayClock [Minute] ) 

{ 

Minute++; 



/* Mode Data Processing */ 
processModes (z, weight, &Cgroup) ; 



/* Calculate delays, shifts, ect */ 
for (k=l; k<=CHANNELS; k++) 

{ 

delay [k] =x [k] [Minute] /Cgroup + /* Time 



delay */ 
Sampling delay 
samples */ 
fraction of 1 



(float) (k-l)*SAMPLE DELAY/1000.; /* 

*/ 

shift [k] = (int) (delay [k] * (float) F_SAMPLE) ; /* # of 

/* 

sample */ 

delta [k] =delay[k] * (float) F SAMPLE- (float) shift [k] ; 

} /* for */ 

} /* if (Need new shift delays) */ 



for (j=l; j<=CHANNELS; j++) 



if (delay [j] >= 0 . 0 ) 



arg= (float) (ORDER+1) /2 . 0+ (1-delta [j] ) ; 
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samples = &inSOUND [j] [i -shift [j] - (ORDER+1) /2] ; 

} /* if */ 

else if (delay [j] <0 .0) 

{ 

arg= (float) (ORDER+1) /2 ,0+delta [j] ; 

samples = &inS0UND[j] [i - shift [j ]- (ORDER+1) /2+1] ; 

} /* else if */ 

INTERPOLATE (indx, samples, ORDER+1, arg, &ans, &err) ; 
buffSOUNDtj] [i-F_SAMPLE] =ans; 
if (ERROR_ESTIMATE==ON) 

MeanSqError [i-F SAMPLE] =MeanSqError [i-F_SAMPLE] +err*err ; 
y /* for */ 

if (ERROR_ESTIMATE==ON) 

MeanSqError [i-F_SAMPLE] =MeanSqError [i-F_SAMPLE] / (float) CHANNELS; 

} /* for */ 

} /* else */ 

/* take fft multiply by weights & sum in frequency domain */ 

df = (float) F_SAMPLE/ (float) FFT_LENGTH; 
if (VERBOSE) printf ("\nTaking FFT\n") ; 
for(i=l; i <=2*FFT_LENGTH; i++) sumXf[i]=0; 

for (j=l; j<=CHANNELS; j++) 

f or (i=l ; i<=FFT LENGTH; i++) 

{ 

if (i <= FFT TIME*F SAMPLE) 

{ 

Xf [2*i-l] =buffSOUND [j] [i] ; 

Xf [2*i] =0; 

} 

else 

{ 

Xf [2*i-l] =0; 

Xf [2*i] =0; 

} /* for i */ 

fourl (Xf , FFT_LENGTH, 1) ; 

k=l; 

for(i=l; i <=FFT_LENGTH-1; i+=2) 

freq=df* (float) (i-l)/2; 
if (freq>=fmin && freq<=fmax) 

{ 

/* positive frequencies */ 

siamXf [i] =sumXf [i] +Xf [i] *weight [j] [k] ; /* real part */ 

sumXf [i+1] =sumXf [i+1] +Xf [i+1] *weight [j] [k] ; /* imag part */ 

/* negative frequencies */ 

SumXf [2*FFT_LENGTH-i] =sumXf [2*FFT_LENGTH- i] +Xf [2*FFT LENGTH-i] *weight [j] 
[k] ; 

sumXf [2*FFT_LENGTH-i+l] =sumXf [2*FFT LENGTH- i+1] +Xf [2*FFT LENGTH- i+1] *wei 
ght [ j ] [k] ; 

k++; 

} /* if (frequency within Bandwidth) */ 

} /* for i */ 
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